Equilibration of Concentrated Hard Sphere Fluids 
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We report a systematic molecular dynamics study of the isochoric equilibration of hard-sphere 
fluids in their metastable regime close to the glass transition. The thermalization process starts with 
the system prepared in a non-equilibrium state with the desired final volume fraction <j) for which 
we can obtain a well-defined non-equilibrium static structure factor So(k; (f>). The evolution of the 
a-relaxation time T a (k) and long-time self-diffusion coefficient Dl as a function of the evolution 
time t w is then monitored for an array of volume fractions. For a given waiting time the plot 
of T a (k; 4>,t w ) as a function of cj> exhibits two regimes corresponding to samples that have fully 
equilibrated within this waiting time (0 < cjy-"' {tw)), and to samples for which equilibration is not 
yet complete {4> > 4>^ c \t w )). The crossover volume fraction 4>^ c \t w ) increases with t w but seems to 
saturate to a value <j>^ = 4>^ c \t w — > oo) « 0.582. We also find that the waiting time t^(<j>) required 
to equilibrate a system grows faster than the corresponding equilibrium relaxation time, t^(cj>) 



0.27 x [ra 5 (k; (j>)] ' , and that both characteristic times increase strongly as (f> approaches 



thus suggesting that the measurement of equilibrium properties at and above ( - a ' is experimentally 
impossible. 

PACS numbers: 05.40.-a, 64.70.pv, 64.70.Q- 



Above a certain size polydispersity, real and simulated 
hard sphere liquids fail to crystalize for volume frac- 
tions 4> beyond the freezing point — 0.494 of the 
monodisperse system As <f> increases the viscosity 

increases enormously, and the metastable liquid eventu- 
ally becomes an amorphous solid. Mode coupling theory 
(MCT) [f| predicts a transition from metastable fluid to 
ideal nonergodic states, characterized by the vanishing of 
the long-time self-diffusion coefficient Dt, and the diver- 
gence of both, the a-relaxation time r Q and the viscos- 
ity r). For the hard-sphere fluid the phenomenology pre- 
dicted by MCT at cf>^ w 0.52 has been essentially con- 
firmed by the experimental observations in hard-sphere 
colloidal suspensions at </>i"p « 0.58 @, 0], although a 
number of intrinsic experimental uncertainties render the 
precise determination of </>i"p a topic of recurrent scien- 
tific discussion 0, Irljioj] . 

The recent work of Brambilla et al. 0, [hJ, however, 
seems to put the very experimental relevance of the diver- 
gent scenario predicted by MCT under severe question- 
ing. By fitting their dynamic light scattering data with 
the asymptotic expression T a (4>) ~ (cj)^ — 0)~ 7 , tradi- 
tionally associated with MCT, these authors determined 
to be ^ = 0.590 ± 0.005. If the ideal MCT picture 
were to be observed in their experiments, the measured 
T a ((j>) should be infinite for <j> > 4>( a K Instead, for the 
volume fraction range (j)^ < <j> < 0.6, they report large 
but finite relaxation times, determined through an ex- 
tremely careful experimental procedure designed to deal 



with artifacts caused, for example, by sample heating or 
sedimentation, which allowed them to accur ately moni- 
tor the equilibration process of their samples 10]. Thus, 
the most immediate interpretation is that these measure- 
ments involve macroscopic states in which the system, 
instead of falling out of equilibrium, remains ergodic and 
enters a new dynamical regime where r Q increases with 
volume fraction according to a different functional form, 
namely, T a (<p) ~ t m expL4(0 o - 0)~ s ]- 

This interpretation, however, rests on the assumption 
that all measurements reporting an apparent stationary 
behavior indeed involve fully equilibrated systems. Re- 



cent molecular dynamics simulations 111, [12[, however, 
suggest that this assumption should not be taken for 
granted without further discussion. For example, ac- 
cording to [11], the relaxation time rtetero of dynamic 
heterogeneities may grow like rtetero ~ T «' 5 as the glass 
transition is approached. Thus, if one has to wait until 
"slow" regions become "fast" regions and vice versa, one 
possibility that cannot be ruled out is that when the equi- 
librium relaxation time T^ q (4>) indeed diverges, the sys- 
tem will require a similarly divergent equilibration time 
t-ivW, thus blurring even the most accurate observation. 
Motivated in part by these considerations, here we in- 
tentionally study the effects on T a ((f>) of the incomplete 
equilibration of concentrated hard-sphere systems close 
to the glass transition by means of systematic computer 
simulations, in which some of the intrinsic uncertainties 
of the experimental samples will be absent. 
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As in 12j, the basic simulation experiment consists 



of monitoring the irreversible evolution of a hard-sphere 
system initially prepared at a non-equilibrium state char- 
acterized by a prescribed volume fraction <p and by a well- 
defined non-equilibrium static structure factor So(k;4>). 
The irreversible evolution to equilibrium is then de- 
scribed in terms of the time-evolving non-equilibrium 
static structure factor S tra (fc; </>) and self intermediate 
scattering function (Sclf-ISF) Fs(k, r, t w ), where t w is the 
evolution (or "waiting") time after the system was pre- 
pared. The naturally expected long-i^, asymptotic limit 
of these properties is, of course, the equilibrium static 
structure factor S eq {k;4>) and self-ISF Fg q (k,r). Our 
interest is to determine the volume fractions for which 
equilibrium is reached within a given waiting time t w . 

We use event-driven molecular dynamics to simulate 
the evolution of N — 1000 particles in a volume V, 
with particle diameters a evenly distributed between 
cr(l — w/2) and cf(l + w/2), with a being the mean di- 
ameter. We consider the case w = 0.3, corresponding to 
a polydispersity s a = w/y/12 = 0.0866. According to 
the results reported in [3] , at this polydispersity the sys- 
tem shows no evidence of crystallization for any volume 
fraction 4> — (ir/6)no~ 3 , where a 3 is the third moment 
of the size distribution and n is the total number density 
n = N/V . All the particles are assumed to have the same 
mass M. The length, mass, and time units employed are, 
respectively, a, M, and a^/M/ksT. 

To produce the initial configurations we used soft- 
particle molecular dynamics to simulate the evolution of 
a set of initially overlapping, randomly placed particles, 
with the correct distribution of diameters, interacting 
through a short-ranged repulsive soft (but increasingly 
harder) interaction, and in the presence of strong dissi- 
pation. For <fi below the random close packing limit, this 
system evolves rapidly into a disordered configuration 
with no overlaps. These non-thermalized hard-sphere 
configurations are then given random velocities generated 
by a Maxwell-Boltzmann distribution, with fcgT set as 
the energy unit. These configurations are then used as 
the starting configurations for the event-driven simula- 
tion of the HS equilibration process. 

The simulations were carried for an array of values 
of <fr between 0.480 and 0.595. For each such volume 
fraction we used waiting times from 1 to 10 5 in pow- 
ers of 10. The sequence of configurations obtained 
is employed to generate the Sclf-ISF Fs(k,T,t w ) = 

(i/^)(Eti ex p[ ik -(^( < » + T )- r »(^))]> and the 

mean squared displacement (MSD) ((Ar(r; t w )) 2 ) = 
(WXEili + r) - r 4 (^)] 2 ), where r^t) is the 
position of the zth particle at time t, r is the corre- 
lation time, and the brackets indicate averaging over 
(at least) 20 independent realizations. Fs(k 7 T;t w ) 
is evaluated at k — 7.1, close to the main peak of 
S eq {k\4>) for all the values of (f> considered. The 




FIG. 1: (a)Self intermediate scattering function Fs(k,r;t w ) 
of a polydisperse hard-sphere system (s = 0.0866) evalu- 
ated at k = 7.1 at volume fraction <f> — 0.575 and poly- 
dispersity s = 0.0866 as a function of the correlation time 
t for waiting times t w = 10°, 10 1 , 10 5 . The inset of (a) 
shows the corresponding So(k; 4>) (black circles) and S eq (k; <j>) 
(red squares). (b) Simulation data of the a-relaxation time 
r a (k; <j>,t w ) as a function of t w at fixed volume fraction. The 
asterisks highlight the points (t^((j>), T^ q (k; <j>)). 



a-relaxation time T a (k; (j),t w ) is defined by the condi- 
tion Fs(k,T a ,t w ) — 1/e, and the long-time self-diffusion 
coefficient Dl by DL{(f)]t w ) = lim T ^oo((Ar(T; t w )) 2 )/6r. 

Let us illustrate the results of this procedure for one 
specific volume fraction, namely, cf> = 0.575. In Fig. QJa) 
we present the simulation results for Fs(k,r:t w ) evalu- 
ated at k = 7.1 as a function of the correlation time r 
for the sequence of waiting times t w — 10°, 10 1 , 10 5 . 
This sequence exhibits the increasing slowing down of the 
dynamics as the system approaches its equilibrium state 
and illustrates the fact that Fs{k,T]t w ) saturates to its 
equilibrium value Fg q (k,r) after a certain equilibration 
waiting time t &q ((f>). For example, from the illustrative 
data in the figure we find that ^((f) = 0.575) « 10 4 . 
A similar evolution and saturation is observed in the 
static structure factor St m (k;4>), which exhibits, as ex- 
pected, a large increase at the first diffraction peak. The 
inset of Fig. [TJa) presents the initial static structure 
factor So(k;4>) = S tnl= o(k; <p) and the final equilibrium 
S eq (k](f). From the data for Fs(k,r;t w ) in this figure 
we can determine the a-relaxation time r a (k; 4>,t w ) as a 
function of t w . The results indicate that the a-relaxation 
time T a (k;4> = 0.575, iu,) saturates approximately to its 
equilibrium value T^ q (k; (j> = 0.575) «2x 10 2 within the 
equilibration waiting time t eq {§ = 0.575) « 10 4 . 

Fig. IHb) plots the dependence of the a-relaxation 
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time r a (k] <fi,t w ) as a function of waiting time t w for 
fixed volume fraction <fi. These plots confirm that 
beyond an equilibration waiting time t^((p), the a- 
relaxation time r a (k;4>) saturates approximately to its 
equilibrium value r^ 9 (fc;0). To emphasize these con- 
cepts we have highlighted the points (t^ q ((f>), T^ q (k; </>)) 
in the figure. In fact, we notice that the highlighted 
points (t^(4>), T^ q (k; </>)) obey the approximate relation 
C(0) ~ °- 27 x [^(fc;^)] 1 ' 43 ^ suggesting that the wait- 
ing time t^((j)) required to equilibrate a system is al- 
ways longer than the corresponding equilibrium relax- 
ation time T^ g (k;<j>), and that both characteristic times 
increase strongly with <f>. 




FIG. 2: Simulation data of the a-relaxation time T a (k; 4>,t w ) 
of a polydisperse hard-sphere system (s = 0.0866) as a func- 
tion of volume fraction at fixed waiting time. The asterisks 
indicate, in this case, the crossover volume fraction <f^ c \t m ) 
at the various waiting times. The circles in the inset display 
the evolution of <j) with waiting time, and the squares 

are the results in Fig. 5(a) of Ref. [12 ] . 



Just like Hermes and Dijkstra [TJ] did with the pres- 
sure (see their Fig. 1), let us display our data for 
T a (k; (f>, t w ) of Fig. [TJb) in a complementary manner, 
namely, as a function of volume fraction for fixed wait- 
ing time t w , and this is done in Fig. [2] The first fea- 
ture to notice in each of the corresponding curves is 
that one can distinguish two regimes in volume fraction, 
namely, the low-0 (equilibrated) regime and the high- 
ly (non-equilibrated) regime, separated by a crossover 
volume fraction 4>^ c \t w ). Focusing, for example, on 
the results corresponding to t w = 10 3 , we notice that 
(j)^(t w = 10 3 ) « 0.57. In Fig. |we have highlighted the 
crossover points ((f>(°'(t w ), r^ q (k; </>)). We observe that 
the resulting crossover volume fraction 4>^ c \t w ) first in- 
creases rather fast with t w , but then slows down consid- 
erably, suggesting a saturation to a value slightly larger 
than 0.58, as indicated in the inset of Fig. [2j which also 
include the results for (j)^ (t w ) determined by Hermes and 
Dijkstra [l2| from the pressure, denoted by rj g in their 
Fig. 5(a). The estimate of the limit <f>^(t w — > oo) will 



hardly be determined by even more powerful simulations, 
and the need of a theoretical framework is required. 

One of the main products of the simulation results just 
presented is the determination of the volume fraction de- 
pendence of the equilibrium a-relaxation time T^ q (k;(j)). 
Clearly, our simulation experiment can determine this 
property only within the window < <fi < <fi^ {t™ ax ), 
where t™ ax is the maximum waiting time achieved in 
the simulation experiment. In our case, t™ ax = 10 s 
yielding (f>^ (t™ ax ) ps 0.58. These results, scaled with 
t q , o (/c;0) = l/k 2 D°, with D° = (see below), 

are plotted in Fig. [3] as solid squares. For > 0.58 the 
^-dependent a-relaxation time T a (k\<fi,t w ) did not sat- 
urate to its equilibrium value within the total duration 
of the present simulation experiment. These results are 
also plotted in Fig. [3] as empty squares, to denote insuf- 
ficient equilibration. Thus, only the data in solid squares 
are meaningful when comparing with the predictions of 
equilibrium theories such as MCT or the more recently 
developed self-consistent generalized Langevin equation 
(SCGLE) theory Q. 

MCT and the SCGLE theory provide similar answers 
regarding the asymptotic divergence of the relaxation 
times. We consider, however, that there is no need to 
appeal to asymptotic expressions, which have a more re- 
stricted range of validity, when one has easy access to 
the full numerical solution of the corresponding theory, 
as we do for the SCGLE theory of colloid dynamics. 
As we have recently discovered jljj] . the latter theory 
also describes the long-time dynamics of atomic systems 
provided the solvent short-time self-diffusion coefficient 
D° is replaced by the kinetic-theory self-diffusion coeffi- 
cient, given by D° = ^/W(p[a^/k^T/M] [H, 0. In 
Fig. [3] we compare the simulated equilibrium data for 
T*(k;<p,t w ) = k 2 D°T a (k;(t),t w ) (<f> < 0.58) with the pre- 
dictions of the SCGLE theory (Ecs. (1), (2), and (5-8) of 
Ref [13j, with k c — 8.48 adjusted to fine-tune the com- 
parison with these equilibrium data). According to this 
fit, T^ q (k-(j>) would diverge at cj>^ sa 0.582. 

There is, of course, no reason to include the non- 
equilibrated data of Fig. [3] in this comparison. As a 
mere fitting exercise, however, we notice that the full set 
including these non-equilibrated data can be fitted by 
the expression r Q (0) = r m exp[^4((/>o — 4>)~ s ], thus find- 
ing A = 0.02, S = 1.921, = 0.21 and <fi = 0.6235 
(dashed line in the figure). Amazingly enough, we find 
that this functional form provides a reasonable fit also 
for the shorter waiting times t w = 10 4 and 10 3 using the 
same values for 5, C, and t^, but with 4>q — 0.631 and 
0.635, respectively, as illustrated in the inset (a) of Fig. 

m 

In order to relate our simulation results with the exper- 
imental observations of Refs. 

dEH, in the inset (b) of 
Fig. [3] we compare the experimental equilibration data of 
Fig. 6 of Ref. for the sample labeled 4> exp = 0.5876, 
with the simulation data corresponding to <p = 0.58 in 
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FIG. 3: Volume fraction dependence of the scaled a- 
relaxation time r* (k; 4>,t w ) = k 2 D°r a (k; 4>,t w ). The solid 
(empty) squares denote simulation data of fully equilibrated 
(insufficiently equilibrated) systems. The solid line are the 
predictions of the SCGLE theory. The dashed line is the fit 
with T a (<p) = Too expL4(0o — 4>)~ S \- The solid circles corre- 
spond to the experimental results of Fig. 13 of Ref. [Io| . 
The inset (a) plots r a (k\ (f>,t w ) as a function of A((f>o — <j>)~ S 
for t w = 10 5 (squares), 10 4 (diamonds) and 10 3 (triangles), 
with the arrows pointing at the corresponding crossover vol- 
ume fraction 4>^ c \t w ). The inset (b) compares our simu- 
lation results (empty squares with line) for t* (k; <f>,t w ) vs. 
t^j = k 2 D°t w for 4> = 0.58 with the experimental data of Fig. 
6 of Ref. OH (empty circles) at <t> exp = 0.5876. 



Fig. HJb) above. The excellent agreement between the 
simulated and the experimental equilibration data sug- 
gests that the difference in the value of 4> and <jf xp could 
be explained by the intrinsic uncertainties discussed in 
Ref. nrj regarding the determination of the volume frac- 
tion of the system. 

Assuming that this is the case, we can directly compare 
our MD simulation results in FigJ3]for T^ q (k; (f>) I T^ q n (k\ 4>) 
with the experimental data in Fig. 13 of Ref. [LOj, 
provided that we assume a constant ratio <f>/(jf xp . The 
dark circles in that figure are precisely those experimen- 
tal data as a function of the experimental volume frac- 
tion reduced by a factor (f>/(f) exp (fitted to yield the value 
0.985), to approximately account for the referred uncer- 
tainties. In addition, as a simple manner to treat hy- 
drodynamic interactions 



17] 



we have to take into ac- 
count that the role of the normalizing parameter D° 
is played, in the experimental data, by the short-time 
self-diffusion coefficient Ds{4>), given approximately by 
D s (<j))/D s (<j) = 0) = (1 - </>)/(l + 1.5<f>) [3. This com- 
parison suggests a completely similar phenomenology, al- 
though it is quite clear only in the case of our simulation 
data that the departure of T a (k;(f),t w — 10 5 ) from the 
equilibrium curve predicted by the SCGLE theory near 
</>( a ) is due to the insufficient equilibration of the system 
within the maximum waiting time t™ ax = 10 5 of our 
simulation experiment. 



The results just presented suggest, however, that any 
simulation aimed at determining the equilibrium value 
of dynamic order parameters such as r a (k; (f>,t w ) and 
Dl(4>] tw) near the dynamic arrest transition is bound to 
be limited by the duration of the simulation experiment, 
represented by the maximum waiting time t w involved. 
This limits the determination of these equilibrium values 
to the window of volume fractions < < 4>( c > (t™ ax ) . 
For 4> > 4>^ (i™ aa: ), the simulation results will be report- 
ing the properties of an insufficiently equilibrated system. 
The results presented here indicate that if we want to en- 
large this window we would have to go to exponentially 
longer waiting times, which is bound sooner or later to 
become a lost battle. There is, of course, no obvious 
reason to believe that a different situation will prevail in 
experimental samples. 

Let us finally notice that the expression T a (<p) — 
Too expL4(</> — c f ) )~ 5 ] gives a reasonble fit for our non- 
equilibrium data, even at early stages in the waiting 
time t w , suggesting that this dynamical regime could be 
consecuence of the lack of equilibration. These means 
that the correct analysis of the data corresponding to in- 
completely equilibrated conditions must be made in the 
framework of a non-equilibrium theory. It is pertinent 
to mention that our original motivation to carry out the 
present simulations was precisely the need to generate 
reliable data incompletely equilibrated systems that will 
serve as a reference to test the recently developed non- 
equilibrium extension of the SCGLE theory. 
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